#include <iostream>
#include<sys/time.h>
#include<immintrin.h>
#include<x86intrin.h>
using namespace std;
const int N=5001;
float m[N][N];

void m_reset(int n){
    for(int i=0;i<n;i++){
        for(int j=0;j<n;j++)
            m[i][j]=0;
        m[i][i]=1.0;
        for(int j=i+1;j<n;j++)
            m[i][j]=rand();
    }
    for(int k=0;k<n;k++)
        for(int i=k+1;i<n;i++)
            for(int j=0;j<n;j++)
                m[i][j]+=m[k][j];
}

void normal(int n) {
    struct timeval h, t;
    cout<<"Normal\n";
    gettimeofday(&h, NULL);
    for(int k = 0; k < n; k++) {
        for(int j = k + 1; j < n; j++) {
            m[k][j] = m[k][j]/m[k][k];
        }
        m[k][k] = 1.0;
        for(int i = k + 1; i < n; i++) {
            for(int j = k + 1; j < n; j++) {
                m[i][j] = m[i][j] - m[i][k]*m[k][j];
            }
            m[i][k] = 0;
        }
    }
    gettimeofday(&t, NULL);
    cout<<"TIME: "<<(1000*(t.tv_sec - h.tv_sec) + 0.001*(t.tv_usec - h.tv_usec))<<" ms"<<endl;
}
//AVX全部并行（不对齐）
void avx(int n) {
    struct timeval h, t;
    cout<<"AVX not align\n";
    gettimeofday(&h, NULL);
    for(int k = 0; k < n; k++) {
        float* va = new float[8];
        *va = *(va+1) = *(va+2) = *(va+3) = *(va+4) = *(va+5) = *(va+6) = *(va+7) = m[k][k];
        __m256 vt = __builtin_ia32_loadups256(va);//将八个单精度浮点数从内存加载到向量寄存器
        delete[] va;
        for(int j = k + 1; j + 7 < n; j += 8) {
            __m256 va = __builtin_ia32_loadups256(&m[k][j]);
            va = _mm256_div_ps(va, vt);//向量化除法运算
            __builtin_ia32_storeups256(&m[k][j], va);//储存到内存
        }
        for(int j = n - n%8; j < n; j++) {
            m[k][j] = m[k][j]/m[k][k];
        }//结尾几个元素串行计算
        m[k][k] = 1.0;
        for(int i = k + 1; i < n; i++) {
            float* aik = new float[8];
            *aik = *(aik+1) = *(aik+2) = *(aik+3) = *(aik+4) = *(aik+5) = *(aik+6) = *(aik+7) = m[i][k];
            __m256 vaik = __builtin_ia32_loadups256(aik);//将八个单精度浮点数从内存加载到向量寄存器
            delete[] aik;
            for(int j = k + 1; j + 7 < n; j += 8) {
                __m256 vaij = __builtin_ia32_loadups256(&m[i][j]);
                __m256 vakj = __builtin_ia32_loadups256(&m[k][j]);
                __m256 vx = __builtin_ia32_mulps256(vaik, vakj);//向量化乘法运算
                vaij = __builtin_ia32_subps256(vaij, vx);//向量化减法运算
                __builtin_ia32_storeups256(&m[i][j], vaij);//存储到内存
            }
            for(int j = n - n%8; j < n; j++) {
                m[i][j] = m[i][j] - m[i][k]*m[k][j];
            }//结尾几个元素串行计算
            m[i][k] = 0;
        }
    }
    gettimeofday(&t, NULL);
    cout<<"TIME: "<<(1000*(t.tv_sec - h.tv_sec) + 0.001*(t.tv_usec - h.tv_usec))<<" ms"<<endl;
}
//AVX全部并行（对齐）
void avxAlgn(int n) {
    struct timeval h, t;
    cout<<"AVX align\n";
    gettimeofday(&h, NULL);
    for(int k = 0; k < n; k++) {
        float* va = new float[8];
        *va = *(va+1) = *(va+2) = *(va+3) = *(va+4) = *(va+5) = *(va+6) = *(va+7) =m[k][k];
        __m256 vt = __builtin_ia32_loadups256(va);//将四个单精度浮点数从内存加载到向量寄存器
        delete[] va;
        for(int j = k + 1; j < k + 8 - k%8; j++) {
            m[k][j] = m[k][j]/m[k][k];
        }//串行计算至对齐
        for(int j = k + 8 - k%8; j + 7 < n; j += 8) {
            __m256 va = __builtin_ia32_loadups256(&m[k][j]);
            va = __builtin_ia32_divps256(va, vt);//向量化除法运算
            __builtin_ia32_storeups256(&m[k][j], va);//储存到内存
        }
        for(int j = n - n%8; j < n; j++) {
            m[k][j] = m[k][j]/m[k][k];
        }//结尾几个元素串行计算
        m[k][k] = 1.0;
        for(int i = k + 1; i < n; i++) {
            float* aik = new float[8];
            *aik = *(aik+1) = *(aik+2) = *(aik+3) = *(aik+4) = *(aik+5) = *(aik+6) = *(aik+7) = m[i][k];
            __m256 vaik = __builtin_ia32_loadups256(aik);//将四个单精度浮点数从内存加载到向量寄存器
            delete[] aik;
            for(int j = k + 1; j < k + 8 - k%8; j++) {
                m[i][j] = m[i][j] - m[i][k]*m[k][j];
            }
            for(int j = k + 8 - k%8; j + 7 < n; j += 8) {
                __m256 vaij = __builtin_ia32_loadups256(&m[i][j]);
                __m256 vakj = __builtin_ia32_loadups256(&m[k][j]);
                __m256 vx = __builtin_ia32_mulps256(vaik, vakj);//向量化乘法运算
                vaij = __builtin_ia32_subps256(vaij, vx);//向量化减法运算
                __builtin_ia32_storeups256(&m[i][j], vaij);//存储到内存
            }
            for(int j = n - n%8; j < n; j++) {
                m[i][j] = m[i][j] - m[i][k]*m[k][j];
            }//结尾几个元素串行计算
            m[i][k] = 0;
        }
    }
    gettimeofday(&t, NULL);
    cout<<"TIME: "<<(1000*(t.tv_sec - h.tv_sec) + 0.001*(t.tv_usec - h.tv_usec))<<" ms"<<endl;
}

int main()
{
    int n=4096;
    m_reset(n);
    normal(n);
    m_reset(n);
    avx(n);
    m_reset(n);
    avxAlgn(n);
    return 0;
}
